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ABSTRACT 

We study the propagation of ultra-high-energy cosmic rays (UHECR) in the 
Galaxy, concentrating on the energy range below the ankle in the spectrum at 4 
EeV. A Monte-Carlo method, based on analytical solutions to the time-dependent 
diffusion problem, is used to account for intermittency by placing sources at ran- 
dom locations. Assuming a source population that scales with baryon mass 
density or star formation (e.g. long GRB), we derive constraints arising from in- 
termittency and the observational limits on the composition and anisotropy. It is 
shown that the composition and anisotropy at 10^* eV are difficult to reproduce 
and require that either a) the particle mean free path is much smaller than a gyro- 
radius, implying the escape time is very long, b) the composition is heavier than 
suggested by recent Auger data, c) the ultrahigh-energy sub-ankle component is 
mostly extragalactic, or d) we are living in a rare lull in the ultrahigh-energy 
cosmic-ray production, and the current ultrahigh-energy cosmic-ray intensity is 
far below the Galactic time average. We therefore recommend a strong observa- 
tional focus on determining the UHECR composition around 10^^ eV. 

Subject headings: Galaxy, cosmic rays, gamma-ray bursts 

1. Introduction 

An open problem in cosmic-ray (CR) astrophysics is at what energy we observe the 
transition from a Galactic to an extragalactic CR origin, and what sorts of sources, in and 
out of our own Galaxy, are responsible for each population. It is widely suspected that 
the abrupt flattening of the spectrum at = 4 x 10^^ eV, the so-called "ankle" of the CR 
spectrum, indicates the Galactic to extragalactic transition, but the abrupt increase in loss 
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rate above the ankle due to pair production with the microwave background may also play a 
role. It is generally recognized that blast waves fr om supernova remnants c annot accelerate 
CR protons beyond energies of about 10^^'^ eV (jLagage fc CesarskyI Il983l ). so they could 
not reach the CR ankle, unless the Galactic m agnetic field is amp li fied b y the blast wave 
by many orders of magnitude. This motivated iLevinson fc Eichlerl ( 1l993l ) to propose that 
UHECR are produced by GRB that occur in our Galaxy. 

This is the third paper of a series, earlier publications of which demonstrated that a) 
long GRB likely have insufficient power to provide UHECR at energies above the observed 
ankle at 4 x 10^^ eV f lEichler et al.ll2010[ ). unless their CR power is much higher than their 
7-ray power, and b) at energies below the ankle the contribution from Galactic GRB should 
on average be stronger by about 3 orders of magnitude than that of extragalactic GRB (as 
the latter would imply that extragalactic space is filled uni formly) , thus further constraining 
models involving a GRB origin of trans-ankle particles (lEichler fc Pohll 1201 ll . henceforth 
referred to as Paper II). 

In this paper we study the time-dependent diffusive transport of UHECR in the Galaxy 
using the method of Monte-Carlo to account for the unknown location and explosion time 
of GRB or other sources with similar population statistics. For simplicity, we shall call the 
sources GRBs throughout the paper, but it is to be understood that the study applies to other 
classes of UHECR sources as well (e.g. 7-ray silent hypernovae). Our approach permits us 
to quantitatively account for intermittency effe cts in the local UHECR spectrum and thus 
goes beyond the s cope of earlier publications (ILevinson fc Eichlerl Il993t IWick et al.l l2004l : 
Calvez et al.ll2010l ). Faced with the observed energy spectrum and low limits on anisotropy 
recently set by Auger and other experiments, we calculate spectra and anisotropy that would 
be expected as a function of model parameters. The ceiling on anisotropy, of order 1 percent 
or less below 10^* eV, is interesting given the rather large mean free paths of such energetic 
cosmic rays. 

We assume the propagation in the Galaxy of cosmic rays at energies 10^^ eV to 10^^ eV 
can be accurately described as isotropic diffusion. The framework of isotropic diffusion 
requires that the particle Larmor radii be considerably smaller than the dimensions over 
which the diffusion is considered. If this were not the case, our conclusion that isotropy is 
hard to understand would be strengthened, as there would be little chance of isotropizing 
UHECR emerging from an anisotropic source distribution. Likewise, the particle mean free 
path, Ainfp, should be much smaller than a few kpc, the typical distance between the solar 
system and a GRB in the Galaxy. The Larmor radius of a Z = 1 particle in a 10 /zG 
field is used as a scale in this paper. It reaches ~ 100 pc at 10^^ eV, and therefore the 
first approximation should hold for UHECRs of any composition below ~ 10^^ eV. The 
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second approximati on requires that Amfp be within a factor of ~ 10 of the Larmor radius. 



Aloisio et al.l (120091 ) have studied the transition from the rectihnear regime to the diffusive 
regime of particle propagation, concluding that diffusion models are a good approximation 
if the source distance is L > 6 Amfp- In the models presented in Section [31 for protons at 
E = 4 EeV Amfp < 500 pc, and smaller by a factor 1/Z for heavy nuclei, whereas the GRB 
rate is so low that typically L > 5 kpc. As we shall demonstrate, reproducing the observed 
low anisotropy in the EeV band requires a mean free path small enough to maintain validity 
of the diffusion approach, at least for GRB. 

Isotropic diffusion requires that the magnetic field is weakly ordered. Measurements of 
Faraday rotation indicate that the strength of the ordered field is about 2 /x G, whereas the 



random field has an amplitude of ~ 5 /iG on scales of 10 - 100 pc (iHanl 120091 . and references 
therein). Note that UHECR sample the turbulent field in the spectral band where they are 
strongest, on the scales of energy injection by supernovae and stellar winds. 

To evaluate the level of systematic uncertainties in our model description, we explore 
various geometric forms of the propagation volume of UHECR in the Galaxy. We find that a 
disk-like geometry, which appears more likely to be accurate than the assumption of spher- 
ical symmetry, renders the observational constraints on anisotropy and composition more 
difficult to meet, and resolution of this problem could lead to interesting or unconventional 
conclusions. 



2. Propagation of UHECRs from Galactic GRB 

2.1. Spherical symmetry 

Treating energy losses as a catastrophic loss term with loss time T, assuming the par- 
ticles are instantaneously released at time to = by a single point source, and assuming a 
large cosmic-ray halo, the problems is spherically symmetric with respect to the GRB, and 
the differential density of UHECR obeys a continuity equation that can be written as 

Here, D = cAmfp/3 is the spatial diffusion coefficient, r is the distance from the GRB, and 
Q{E) is the differential production rate of UHECR in a GRB. We assume that we are located 
outside of the dissipation region of the GRB, implying that the plasma no longer streams with 
significant velocity. If we were, the Compton-Getting anisotropy would be prohibitively large. 
Once the GRB outflow turns non-relativistic, the further dynamical evolution is not unlike 
that of a supernova remnant, and we can expect the GRB remnant to be of order 100 pc in 
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size, small enough to render ii highly unlikely that we reside inside. The high-energy particles 
escaping from a GRB ca n initially have a rather directional distribution wh ich leads to well- 
known instabilities (e.g. iMedvedev fc LoeblllQQQt iPohl &: Schlickeiserll2000| ) that provide for 



a rapid isotropization of the particles. We assume that this isotropization is achieved on a 
scale small compared with the average distance to GRBs, and hence we can approximate it 
as being instantaneous. 

Note that equation [T] contains the particle energy, E, only parametrically, it is not a 
variable. The catastrophic loss term permits using Dirichlet boundary conditions N{r = 
oo) = 0. 



Using standard methods (lKardashevlll962l : iGinzburg fc Syrovatskiilll964j ) one derives a 
general solution to the transport equation, 

^2 



N{r, t, E) = exp 
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■Q{E) exp 



4:Dt 



(2) 
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If more than one GRB contribute to the local UHECR flux at any time, their individual 
contributions must be calculated using equation [2] and then summed. 

The anisotro py in the UHECR intensity mainly arises from the diffusive flux (e.g. 
Schlickeiserl Il989[ ) and is of the order 
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where A^tot is the total UHECR spectrum after summing the contributions of all relevant 
GRB. If there was only one GRB, then the anisotropy would be 

' (4) 
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and therefore, setting for example r = r0, a given anisotropy at the solar circle, 5, requires 
that the age of the UHECR population, t, obey 
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At that time, the Galactic UHECR fill a sphere of radius 
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(6) 



100 pc 

To be noted from the approximation represented by equations H] and [6] is that the UHECR 
flux at the solar circle depends on the mean free path, but the anisotropy does not (e.g. 
Wick et al.l[200i ). 
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The above conclusion is highly significant. It means that at increasingly high energies 
the anisotropy may be reduced, given a sufficiently long time since the last point-source 
injection, to within the tight observational constraints even at energies high enough that the 
mean free path is a significant fraction of the geometric scale of the Galaxy, as observations 
indicate. We therefore investigate how the solution generalizes to arbitrary radial dependence 
of the diffusion coefficient. 

Consider, then, equation H] with a radially dependent diffusion coefficient D = DqT^ , 
with j > 0. 

^ + ^ - fr^ ^1 = QiE) 6(t) ' (7) 

Note that we choose a finite radius of injection. To, in a way that preserves spherical symme- 
try, so that the diffusion coefficient at the injection site is finite. For j < 2, a finite boundary 
condition at the origin and a vanishing boundary condition at r = oo, the solution at large 
times (or r ^ Tq) approaches 

N{r, t, ~ exp — Q{E) exp ' 



By comparison with Equation [21 the case of constant diffusion coefficient, the solution falls 
off more quickly after reaching a peak ffux at = 3 (2 — j) Dq t, whereas the anisotropy, 
when estimated as in Equation [HI is 

(which differs from the expression in Equation [H by the numerical factor of order unity 
(1 — |)). That the anisotropy has similar time dependence, while the decline in intensity has 
steeper time dependence has further significance: It means that to reduce the anisotropy with 
time at the observer's location, a larger reduction in intensity is required for positive j than 
for J = (the case of spatially homogeneous diffusion). If we were to attempt accommodating 
the low observed UHECR anisotropy with an anisotropic source distribution by invoking an 
intermittent source of UHECR followed by a period of zero source activity, we would have 
to pay a larger cost in intensity reduction. The intensity at Earth at the time of maximum 
intensity, we would conclude, would be higher than the current intensity by an even larger 
factor than in the case of homogeneous diffusion. 



When j > 2 the exponential factor does not vanish at r = oo, and the correct Green's 
function at large time is 
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(10) 

In this case, the scale height of N at radius r, (din A^/rfr)"^^, is of order r, even at large t, 
and the anisotropy does not become arbitrarily small in the limit of large t. If the magnetic 
field of the Galaxy is dragged out by a Galactic wind, it is quite reasonable to suppose that 
the magnetic field decreases as and that the gyroradius and mean free path increase 

at least as fast as r^. It is therefore not clear whether the small observational upper limit 
on anisotropy can be accommodated with a realistic magnetic field geometry, no matter 
how intermittent the sources. Similar remarks would apply if convection were significant as 
compared to diffusion. Moreover, we have not yet taken into account the fact that a single 
intermittent source of cosmic rays would be localized in angle as well as radius, which would 
further increase the anisotropy. 

Below we reconsider the anisotropy question assuming a free-escape boundary up to 
several kpc above the Galactic disk following the consideration of several other issues. The 
free-escape boundary models the effect of a very strong dependence of D with r (e.g. due to a 
weaker magnetic field) and/or the outward convective effects of a Galactic wind or of Galactic 
chimneys beyond some region. We show that the same principle applies as demonstrated 
above in the spherically-symmetric case with spatially varying coefficient: that the leakage 
out the faces of the disk greatly reduces the intensity of the UHECR as they isotropize 
relative to the observer; in fact, while the anisotropy declines algebraically with time, as in 
the previous examples, the intensity declines exponentially, so the isotropization carries a 
large cost in intensity reduction. This poses constraints on any theory of their origin from 
sources that are distributed anisotropically in the observer's sky. 



2.1.1. Heavier nuclei 

Heavier nuclei have a smaller rigidity than protons at the same total energy, R cc E/Z. 
The mean free path of an ultra-high-energy particle should only depend on the rigidity, and 
in the absence of energy losses a nucleus of charge Z and energy Ez should behave like a 
proton of energy E = Ez/Z. Thus equation [2] also describes the distribution of heavy nuclei 
in the Galaxy, provided the appropriate scaling is applied to the energy and the source rate. 
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2.2. Diffusion in disk geometry 



Cosmic-ray studies at lower energies indicate that a few kpc out into the halo diffusive 
transport become much faster, and particles efficiently escape from the Galaxy. For a cosmic- 
ray source near, e.g., the Galactic Center we therefore expect deviations from spherical 
symmetry, and a thin-disk geometry appears better suited to describe the transport of cosmic 
rays. The cosmic-ray continuity equation ([1]) then changes to 

where the halo boundary can be incorporated through boundary conditions N(z = ±if) = 0. 
Integrating over z yields 



dM M d (jjdM 
dt T dx \ dx 
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Instead of a computationally expensive full solution (jBiisching et al.l l2005l ) . we use a 
steady-state solution to estimate the relation between M and the mid-plane cosmic-ray 
density, Nq, as well as turn the diffusive flux at the halo boundaries into a simple catastrophic 
loss term. In the energy band of interest, escape is the dominant loss process of cosmic rays 
in the Galaxy, and in the steady state we expect their density to follow 
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The diffusive escape can thus be described using a loss time 

r2 / X 2 
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(14) 



(15) 



(16) 



which is far shorter than th e energy-loss time, T, which is ^10^ years for photo-secondary 
production in the Galaxy (IBerezinskii fc Grigoreyal Il988l ) and ~ 10^ years for secondary 
production in p-ISM collisions ( jHuang et al.l 120071 . with average density tt-ism < 0.2 atoms 
/cm^), and therefore energy losses can be ignored henceforth. Ignoring variations in the 
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Fig. 1. — CR density as function of r (p) for isotropic diffusion (solid line) and disk geometry 
(dashed). The parameters are t = 3 Myr, H = 5 kpc, and A = 0.1 kpc. 



diffusion coefficient within the Galactic plane, the problem only depends on the in-plane 
distance between source (GRB) and observer, p, and can be recast as 2-D diffusion equation 
for the mid-plane cosmic-ray density around a point source. 



&t Tesc p dp 



dp 



27C pH 



(17) 



with solution 



iVo(p,t,i.)^exp(-A^^^exp 



The anisotropy in the case of a single GRB is the same as in equation |U using p in lieu 
of r. The ratio between the cosmic-ray density calculated for a disk geometry and that for 
spherical symmetry is 

' ' t 



R 



t 

277 exp 



(19) 
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which, as expected, becomes small if t ^ Tesc- A comparison of the radial dependence of the 
particle density for spherical symmetry as opposed to disk geometry is shown in Figure [1] 
To be noted is that the radial distribution is unchanged, the only difference is seen in the 
normalization. 



2.3. Discreteness Anisotropy 

The above discussion considered anisotropy from a single source in conditions of high 
intermittency. When many sources contribute, the anisotropy is reduced by a factor of order 
the square root of the number of contributing sources. We can formally define intermittency 
for the case of a disk geometry with absorbing boundaries as follows: The CR escape time 
from the disk, Tesc, if the CR are produced in the Galactic plane, is Xesc = H'^/{2D) = 
3if^/(2Amfpc). Here D, the CR diffusion coefficient, is written as Amfpc/3, where Amfp is the 
CR mean free path. The intermittency factor / can be defined to be 1/Mcr where Mcr is 
the number of CR episodes per diffusion area over one escape time. Because the CR escape 
time is the time required to diffuse the height of the disk, H, the diffusion area over one 
escape time is tcH'^, so Mcr = qr^scH'^ /R^-, where R is the radius of the Galaxy and q is the 
rate of CR production events per the entire Galaxy. (For supernovae, gsN is thought to be 
about 1/30 yr, and for GRB, gcRB is thought to be somewhere between 30,000 years and 
3 million years.) Sources much further away than H are attenuated by escape and do not 
contribute significantly. Finally we write Mcr as 

M - ^^^^ (20) 

2 i?2 Amfp c 

and, assuming that the average distance to the sources is of order H and using equation (3), 
we write the anisotropy in the case of M contributing sources as 



hi ^ = ( ^V^' ^ J^. (21) 

H./M^ \ H J ny 3qH ^ ' 

The exact anisotropy of course depends on time and numerical factors of order unity. It 
may be increased by an inhomogeneous spatial distribution of CR sources in the Galaxy, 
and therefore in subsequent sections we shall consider the problem with Monte-Carlo gen- 
eration of CR events. It can be seen however, that 6m, at EeV energies, can be consider- 
able. For a Galactic magnetic field of 10 /iG, the gyroradius of a 3-EeV protons is 0.3 kpc, 
and thus the factor [Amfp/-f^]^''^ is more than 1 percent even if A^fp is but a single gyrora- 
dius. The factor R/H is considerably larger than unity except for a rather thick disk, and 
V 3^17 least of order unity for GRB, and about 10~^'^ for SN. Unless H/R is chosen 
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to be of order unity, the diffusion is sub-Bohmian, or the trans-EeV CR are heavy nuclei, 
the anisotropy exceeds the experimental l imit of 1 percent set by the AUGER experiment 



( iPierre AUGER Collaboration et al.ll201l[ ). The inverse cubic dependence of 5 on if con- 
strains H to be larger than 3 kpc if SN are the sources of UHE Galactic CR, and, if GRB 
are the sources, appears to require them to be intermittent. 

The qualitative point is that, if the sources are distributed in the Galactic plane, then 
" Olbers paradox" does not apply and close sources contribute as much as distant ones. This 
allows the possibility of discrete anisotropy at ultrahigh energies. 



2.4. Intermittency 

2.4- 1- Spectra 



Generally, GRBs in the Galaxy are expected every million years or so, the exact rate 
depending on the beaming fraction and the detailed scaling of lon g GRB with star formation 
and metallicity (For a detailed review see iGehrels et al.l l2009l ). Therefore, only a small 
number of GRB can contribute to the particle flux at the solar circle, and their relative 
contribution depends on the location and explosion time of the GRB. Variations in the local 
particle flux must be expected, and neither the particle spectrum from an individual GRB 
nor the spectrum calculated for a homogeneous source distribution are good proxies. To 
fully account for discreteness of GRBs in space and time, we can use the method of Monte- 
Carlo to randomly place GRBs in the Galaxy with given spatial probability distribution in 
galactocentric radius. 



P{rGc) 



2rGc 



exp 



GC 



with scale tq 



5 kpc, and with given GRB rate, e.g. 

Pit) = 0.5 MjT-\ 



(22) 



(23) 



The assumed galactocentric distribution is an approximation to the distributions of baryons 
and star formation in the Galaxy, which are reasonable proxies for the distribution of super- 
novae. The distribution of lo ng GRB may be skewed toward the out er Galaxy by metallicity- 
related selection effects (e.g. iLevesque et al.ll2010l : iLevesqud 1201 ll ). In this paper we keep 
fixed the galactocentric distribution of UHECR sources. In a forthcoming study we shall 
explore the impact of sources in the outer Galaxy, e.g. short GRB. 

Each GRB is assumed to inject particles over 10^ yr with the same total energy content 
and spectrum 

Q{E) = QoE-' (24) 
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The actual choice of spectral index is irrelevant for the calculation, because it propagates 
right into the final expressions for the particle spectra (cf. equations |2] and [T8|) . Likewise, it 
is sufficient to calculate the variations in the local flux for protons, because nuclei of charge 
Z and energy Ez behave like a proton of energy E = Ez/Z (cf. section [2.1.ip . The halo 
size is again H = 5 kpc. 

Below we show results for two different energy dependences of the diffusion coefficient, 
or mean free path. 

Case A : A^fp = lOr^ (25) 

Case B : A^fp = 10 E^°-^ (26) 

where the Larmor radius is calculated for a proton in a 10-yuG magnetic field, and En = 
E/{10^'^ eV). Case B thus uses a diffusion coefficient that increases with the square root of 
the energy, similar to that conventionally used for Galactic cosmic rays below the knee. Note 
that for case A the diffusion approximation breaks down at approximately 10^^ eV (10^^ eV 
for case B). 

For both case A and case B we have calculated spectra for 10^ random sets of GRBs. As 
we demonstrated before, energy losses through inelastic collisions can be ignored. To reach a 
quasi-steady state, we need to sample the contributions of all GRBs for a time period much 
longer than the escape time. We chose to integrate over 6 Gyr in time, which is long enough 
at all energies and for all choices of Amfp used in this study. Overall, the m ethod is the same 



as that used to model the transport of cosmic-ray electrons in the Galaxy (jPohl &: Esposito 



19981 : iPohl et al.ll2003l ). Results are shown in Figure [2] for case A, and in Figure [3] for case B, 
both with injection index s = 2. The shaded areas in the Figures indicate the central 68% 
containment region for the particle flux, meaning 16% of the simulated spectra fall below 
the colored range, and another 16% are higher. 

To be noted from the flgures are 



• The energy dependence of the diffusion coefficient influences the particle spectrum 
observed at Earth. Features in the observed spectrum could thus arise from changes 
in the energy dependence, e.g. from shallow at lower energies to Bohmian at higher 
energies, without requiring any structure in the source spectrum (see also Calvez et al. 
2010). 

• The solution for disk geometry is typically a factor 2-3 lower than that for spherical ge- 
ometry, reflecting the fact that the free-escape boundary in the halo introduces particle 
losses. 
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Fig. 2. — Proton spectra at the solar circle expected for a diffusion coefficient with Bohmian 
energy scaling (case A). For spherical geometry, the solid line shows the time-averaged par- 
ticle spectrum, whereas the dotted line denotes one randomly selected case. The blue band 
indicates the central 68% containment region for the particle flux at the given energy. The 
dashed line, dashed-dotted line, and the red band (partially obscured) are the same for disk 
geometry. 
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Fig. 3. — Proton spectra at the solar circle expected for a diffusion coefficient scaling with ^J~E 
(case B). For spherical geometry, the solid line shows the time-averaged particle spectrum, 
whereas the dotted line is one randomly selected case. The blue band indicates the central 
68% containment region for the particle flux at the given energy. The dashed line, dashed- 
dotted line, and the red band (partially obscured) are the same for disk geometry. Note that 
the y-axis is now E"^'^ I. 
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Fig. 4. — Proton spectra at the solar circle expected for a diffusion coefficient scaling with 
a/E (case B) in disk geometry. The parameters are as in Figure [3l except that we have 
increased the GRB rate to 5 Myr~^. The thin solid lines show 31 randomly selected spectra, 
whereas the time average is now given by the thick white line. 
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Intermittency is strong for a GRB rate below 1 per Myr, in particular for the more 
realistic disk geometry. In essence, the local UHECR spectrum from Galactic GRBs 
is unpredictable if the scattering mean free path exceeds about 100 pc, which for the 
parameters used here is the case above 10^'' eV for protons, and above 3 ■ 10^^ eV for 



iron. Model fits of single-source spectra can thus be very misleading (cf. IWick et al. 



Particularly interesting is the finding that the time-averaged spectrum can be outside 
the 68% containment region, which arises from the very rare GRB that are both recent 
and near. It indicates that it is unlikely that we observe the time-averaged spectrum. 
The median of the distribution shows a spectral cut off, which is therefore the most 
likely outcome, whereas the time-averaged spectrum does not and thus is not useful in 
a comparison with data. 

The actually expected spectra display bumps unr elated to both source and propa gation 
physics, some of which may indeed be observed fl Arteaga- Velazquez et al. 2010). The 



absence of very large bumps in the observed UHECR spectra suggests that either the 
mean free path for scattering is smaller than assumed here, or the rate of cosmic-ray 
producing GRBs in the Galaxy exceeds 1 per Myr, at which the amplitude of such 
bumps becomes smaller (see Figure H]). Generally, careful accounting of the statistical 
fluc tuations is mandat ory for proper estimating the local UHECR spectrum from GRBs 



fcf. ICalvez et al.l 120101 ) 



2.4-2. Anisotropy 

If more than one GRB contributes to the local UHECR spectrum, their anisotropy 
signal will generally not line up and can thus not be simply averaged. If all sources reside in 
the Galactic plane, the entire problem is essentially two-dimensional. We use an Eddington 
approximation for the UHECR intensity from a single GRB, 

I = Iq + h cos9 = Io{l + 5 cos9), (27) 

where 6 is the angle relative to the line of sight from the GRB to us. The GRB is located 
at Galactocentric coordinates (rQccp), which translates to distance, r, and line-of-sight 
direction relative to the anticenter, ip, as 

r2 = 64 + r^c - 16 tqc cos (p (28) 

and 

8 - rcc cos rcc sin 
cos-^ = smtp = (29) 
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where the distance unit kpc is used. The local anisotropy can be expressed through projection 
in anticenter direction, e^^, and in a direction perpendicular to it, e^. 

Ix = 5 Iq cos ip ly = 6 Iq sin ip (30) 

After summing over all sources, 

-^tot = -^0 -^tot.x = Ix -^tot,y = ly (31) 

the observed anisotropy is 



and the direction is given by 



Stot = ^ (32) 



a = arctan ( ) (33) 

'tot.x 



Figure |5] shows the anisotropy signal for the parameters of case B (see Equation [26!) . For 
Bohm-like scaling of the mean free path (case A), the anisotropy is the same at 10^^ eV but 
has a steeper energy dependence. The current experimental limits on 5tot set by the Auger 
Observatory are 6 t r,t <0.01 in the [0.4,4] EeV range and Stct < 0.05 in the [4,40] EeV range 
JArmengaudlboOsI : bhimer et aPboOol : iBonino et allbood : IPierre AUGER, Collaboration et al" 



20111 ). To be noted from the figure is that the anisotropy for protons above 10^^ eV is high 
compared with the Auger limits, in particular for the more realistic disk geometry. Note 
that the average anisotropy (dashed and solid lines) is independent of the GRB rate. A sub- 
stantial contribution of heavy nuclei, a small mean free path, or extragalactic UHECR would 
be needed to satisfy the anisotropy limit. Figure |6] displays the variation in the anisotropy 
direction, suggesting that one may underestimate the anisotropy if one only considers its 
projection on the anticenter direction, fairly independent of the choice of parameters. 



3. Combined analysis 



We now try to construct a model that reproduces the spectrum of cosmic rays between 
4 ■ 10^® eV and 4 • 10^^ eV together with the a nisotropy limits and the composition. We 



use d ata of the Kascade- Grande collaboration (lApel et al.l l2009l : lArteaga- Velazquez et al. 



20101). ranging frorn 2 PeV to 0.5 EeV, HiRes mo nocular data between 0.2 EeV and 3 EeV 
(lAbbasi et al.ll2005h . and the Auger collaboration ( jPierre AUGER Collaboration et al.ll2010l ) 
above 1 EeV. 



At 10^^ eV the anisotropy is low, 6 < 0.01 (the 99% upper limit is 0.014) flPierre AUGER Col 



201ll ). and the composition is light, though not necessarily dominated by protons flAbraham et al 



aboration 
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0.01 0.10 1.00 10.00 

Energy in lO'^eV 



Fig. 5. — The expected anisotropy for the same parameters as in Figure [31 For spherical 
geometry, the sohd hne shows the average anisotropy, whereas the dotted hne is one randomly 
selected case. The blue band indicates the 68% containment region for the anisotropy at the 
given energy. The dashed line, dashed-dotted line, and the red band (partially obscured) are 
the same for a disk geometry. 
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0.01 0.10 1.00 10.00 
Energy in lO^^eV 



Fig. 6. — The distribution of anisotropy directions relative to the anticenter for disk geome- 
try. The red and blue lines mark the boundaries of the 68% containment range for shallow 
(a=0.5) and Bohmian (a=l) energy scaling of the mean free path, respectively, both for a 
low GRB rate of 0.5 Myr~^. For comparison, the black line is for a=0.5, but high GRB rate 
5 Myr"\ 
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Fig. 7. — Example of model spectra for cosmic-ray protons, helium, and carbon nuclei, 
including the 68% variation range in the case of protons. The solid line denotes the time- 
averaged spectrum of protons, the dashed line displays the same for helium, and the dash- 
dotted line is for carbon. The mean free path follows equation |3ll and the corresponding 
escape timescales are indicated at the top axis. Also shown are the spectra measured with 
KAS CADE- Grande below 0.5 EeV, Auger above 1 EeV, and HiRes in between, together 
with labels indicating the composition. The offset between spectra from different experi- 
ments is likely due to errors in the absolute energy scale. For comparison, we also display 
correspondingly labelled simple mock spectra of protons, helium, and C-Fe-group particles 
possibly produced in SNR. The thick triple-dot-dashed line indicates the total of protons 
from GRB and all particles from SNR. 
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Fig. 8. — Example of model spectra for cosmic-ray protons in the presence of an extragalactic 
proton component. The solid line denotes the median spectrum of protons, and the dash- 
dotted line displays the extragalactic proton spectrum. Otherwise, everything is as in Figure 

m 
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Fig. 9. — Top panel: 200-Myr lightcurve of the local proton flux (cf. Figure [7]) in arbitrary 
units at 3 different escape times relative to the inverse source rate of 1 Myr: The dotted 
(blue) line is for t^sc = 20 Myr , the dashed (cyan) line for t^sc = 4.5 Myr, and the solid 
(red) line for Tesc = 0.5 Myr. The two black bars on the right of the panels, extended as 
thin dashed lines, indicate the median and the average of the light curves at Tesc = 5 ■ 10^ 
yr, where the lower bar always indicates the median. Middle panel: Temporal variation of 
the anisotropy at the same escape times. For the diffusion parameters chosen in the text, 
where the mean free path is chosen conservatively as the gyroradius for the two highest 
energies, the dotted (blue) line is for 0.024 EeV, the dashed (cyan) line for 0.24 EeV, and 
the solid(red) line for 2.4 EeV. The red bar on the right indicates the 99% upper limit set 
by Auger. Bottom panel: Temporal variation of the anisotropy for the model including an 



- 22 - 



2010l ). Figure [5] suggests that the mean free path must be relatively small to accommodate 
the anisotropy limits for light nuclei such as helium. Cosmic-ray propagation models can 
reproduce the composition and particle spectra in the GeV-TeV band for a shallow energy 
dependence of the diffusion coefficient (D oc E^-^^), if stochastic reacceleration is included 
( IStrong fc Moskalenkdll998t iMaurin et al.ll200ll : iJones et al.ll200ll ) . If we were to extrapolate 
their results, the proton mean free path at 10^^ eV would be ~ 500 pc, i.e. larger than what 
we have assumed above. Obviously, we cannot permit the extrapolation to hold, if Galactic 
GRB are to account for the observed UHECRs in the 10^* eV energy band. The mean free 
path must be smaller than the extrapolation, but how much smaller, in particular how small 
compared with the particle Larmor radius? 

Figure [7] shows spectra for a possible model configuration, where for simplicity we display 
only spectra for protons, helium, and carbon as proxies for light and heavy nuclei, respec- 
tively. The GRB rate is set to P(t) = 1 Myr^^ and the source spectral index is s = 2.1. The 
mean free path transitions from a shallow energy dependence to Bohmian scaling (oc E) as 
the particle energy increases. 



A 



mfp 



Ao E' 



0.3 



1 + 



E 



60 PeV 



0.7 



(34) 



where Aq is chosen so a proton has a mean free path that is a certain multiple (unity in 
Figure [7j) of 11 pc at 10^^ eV, its Larmor radius in a 10-/iG magnetic field. The top axis 
of Figure [7] indicates the escape time at a few energies. The spectrum below 10^^ eV is far 
below observed values to accommodate other Galactic sources of cosmic rays, such as SNR 
or PWN, for which we plot simpl e pow er laws with exponential cut-off at Z ■ (3 PeV) to 
guide the eye, based on lAhn et al.l (120101 ). For that we assume the averages < A >= 25 and 
< Z >= 12.5 for the C-Fe group and the following power-law indices: s = 2.66 for protons, 
s = 2.58 for helium, and s = 2.56 for C-Fe. It is apparent that the energy scales of HiRes 
and Auger are discrepant. To fully fit the data we need to increase the normalization of 
the GRB-particle spectrum and either disregard the Auger data or, alternatively, distrust 
HiRes and require a fluctuation in the local spectrum that makes the GRB-produced protons 
somewhat softer than the average spectrum plotted in Figure [71 Such a deviation would be 
common, because the fluctuation amplitude at escape times below 1 Myr is large for the 
GRB rate used here, one per Galaxy per Myr. 

The required spectral shape above 1 Ee V depends on the prope rties of the extragalactic 
component that dominates above the ankle. iBerezinsky et al.l (120061 ) have demonstrated that 
the effect of the pair-production threshold on an extragalactic component composed solely 
of protons can match the ankle fairly well. In Figure M we add such a component using a 
simple approximation of the pair-production loss rate {E oc EQ{E — -Ethr.)? assuming no 
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evolution, i.e. m = in (IBerezinsky et al.ll2006[ ). and integrating out to a redshift Zmax = 5) 



The evolution and redshift c ut-off can in princ iple be constrained using the neutrino and 



gamma-ray backgrounds (e.g. IWang et al.ll201ll ). but such an endeavor is beyond the scope 



of this paper. In calculating the anisotropy (cf. Figure M>, we assume the extragalactic 
component to be perfectly isotropic. The motion of the Sun and the Galaxy relative to the 
CMB frame suggests the pr esence of a Compton-Getting anisotropy of amphtude ~ 0.5% 



(IKachelriefi fc Serpicdl2006l ). but it is unclear to what degree scattering in the Galaxy can 
reduce the observed anisotropy. Assuming no anisotropy in the extragalactic component 
therefore constitutes a best-case scenario (in the sense of diluting the anisotropy) that is 
likely too optimistic. Note that we slightly rescaled the component from galactic GRB and 
plot the median instead of the time-average. 

The average spectrum mainly depends on the absolute value of the mean free path at 
high energies. We can calculate the cosmic-ray source power required to sustain the observed 
flux of UHECRs at 10^^ eV as function of the mean free path. For the parameters used here, 
and an injection spectrum oc E'"^'^ extending from the GeV band to the highest energies, 
fitting the observed flux of UHECRs at 10^^ eV requires the source power 



CR 



P{t) [ dE EQoE-' {W^' erg/ s) (35) 

Jl GeV 



The source power in the energy interval [10 , 10 ] eV alone is 



A 



EeV 



™fP p; in35 



(3.5- 10^^ erg/s), (36) 



which is about 10~^ of the total CR luminosity of the Galaxy. By way of comparison, the 
source intensity of alternative sources to SNR, if they are to supply all Galactic CR at 
energies above the Lagage-Cesarsky limit of E'lc = 10^^'^ eV for SNR, must be the Galactic 
luminosity of CR above these energies, i.e. 

P ~ in(-o-'^)(5-5) r resc(GeV) nn^e-s pro-/«^ ^(GeV) 

PsNR^lO (105.5 GeV) ^esc(105-^GeV)- ^^^^ 

Assuming that Tesc{GeV)/Tesc{lO^-^GeV) > 30, this is about lO^^ erg/s above lO^-^ GeV. On 
the other hand, supernovae may supply the observed CR all the way to 10^*^ eV if those at 
the highest energies are heavy nuclei. 

What Figures [7] and [S] serve to demonstrate is that a GRB produced cosmic-ray compo- 
nent can fill the gap between the contributions of ordinary galactic cosmic-ray sources and 
the super-ankle region likely supplied by extragalactic sources. The discreteness fluctuations 
are typically larger than experimental uncertainties. But what of the anisotropy and its 
variation? 
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Fig. 10. — Top panel: 200-Myr lightcurve of the local proton flux at the 3 different escape 
times denoted in Figure [9l assuming the protons originate in SNR with the parameters 
described in the text. Bottom panel: Temporal variation of the anisotropy at the same 
energies. 
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The lightcurve shown in the top panel of Figure [9] shows that besides the increase 
in fluctuation amplitude, the fluctuation timescale decreases with particle energy (or more 
precisely escape probability). For the parameters used here, the escape time at 2.4 EeV 
is about 5 ■ 10^ yr. The flux variations can be compared with the temporal behaviour of 
the anisotropy, shown in the middle panel of Figure M Already at 0.24 EeV (the dashed, 
cyan line) the anisotropy hovers around 1%. A clear correlation between flux and anisotropy 
is not discernible at higher energies, although at 0.024 EeV some structures in flux and 
anisotropy coincide. However, other spikes in flux or anisotropy observed at 0.024 EeV 
have no counterpart at higher energies, and vice versa. To be noted from the figures is that 
significant excursions to low anisotropy are rare. The minimum anisotropy is 0.0045, whereas 
the minimum flux is 8 • 10~^. 

Figure M also indicates the median and the time average of the flux and anisotropy at 
2.4 EeV. The average is what we would expect to observe if the number of UHECR sources 
were much larger. In contrast, the median is a typical result in the case of intermittent 
sources. To be noted from the figure is that both the average and the median at 2.4 EeV are 
considerably in excess of the Auger limits. Nevertheless, the anisotropy falls below the 99% 



uppe r limit of 2.3% (For the energy band between 2 and 4 EeV lPierre AUGER Collaboration et al 



201ll ) during about 25% of the time. The bottom panel of Figure |9] shows the anisotropy 
variation in the energy band [1, 2] EeV for the case including an extragalactic proton com- 
ponent as in Figure [H which we assume to be perfectly isotropic. The average anisotropy is 
slightly reduced, but we note that the anisotropy becomes very small during lulls, essentially 
because during these times the extragalactic component dominates the spectrum. Still, the 
anisotropy falls below the 99% upper limit of 1.4% only during 40% of the time. 

During the long lulls the intensity dips as low as three or four orders of magnitude 
below the time average! This possibility would suggests an anthropic scenario whereby we 
live in highly unusual times, and that the present cosmic-ray intensity is, perhaps, far below 
usual conditions, which for some (possibly obscure) reason would be hostile to intelligent life 
and/or advanced civilization. Note that the extreme versions of this scenario, in which the 
present UHECR luminosity of the Galaxy is ~ 10~^ below the average (and possibly an even 
larger factor below the most recent maximum) , would require an average (or peak) luminosity 
above 100 GeV or so that compares with (or exceeds) the present day value, and hence an 
extremely flat Galactic source spectrum. This could possibly be tested or constrained by 
abundance measurements of terrestrial cosmogenic nuclei (TCN). It should be remembered 
that the hypothesis of sharp time dependence applies only to the CR that escape rapidly 
enough, and the anisotropy bounds motivate this only at the highest energies. Thus, the 
TCN production rate as a whole is not necessarily strongly affected by the hypothesis, so 
TCN anomalies are not firmly predicted by it. On the other hand, the production of ionizing 
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particles at sea level per unit primary energy is an increasing function of the latter (assuming 
shower maximum occurs above sea level). Thus, for a perfectly flat spectrum (equal energy 
flux per decade of primary energy), most of the ionizing secondary flux at sea level, which 
might possibly be of astrobiological relevance, would come from CR at ultrahigh energies. 
Further discussion of an anthropic scenario is beyond the scope of this paper, but will be 
the subject of future work. It is mentioned here to justify consideration of a scenario that 
would otherwise seem a priori improbable, but accounts for the remarkably low limits on 
anisotropy at ultrahigh energy. 

Also to be noted is that the high anisotropy is unavoidable in scenarios involving UHECR 
acceleration in SNRs, because intermittency is not an issue and there is little uncertainty 
concerning the spatial distribution of SNRs in the Galaxy. Lightcurves for an SNR origin 
are shown in Figure [TUl for which we assumed a supernova rate of one per 50 years. The 
rare spikes in intensity are caused by very nearby supernovae, for which their spatial and 
temporal evolution becomes an issue. Since we maintain the duration of the particle injection 
per source, 10~^ Myr, the spikes are likely clipped by temporal smoothing. 

The anisotropy can in principle be reduced by choosing a smaller mean free path. Figure 
[TT] shows the anisotropy for protons and a mean free path following equation [3H but with 
Ao chosen so that at high energy Amfp = 0.2 r/,. The Auger limits, 5 < 0.01 at 1 EeV, can 
be satisfied with protons. The margin is not large, though, and a few more years of Auger 
operation should permit detecting an anisotropy signal, even for an unusually small mean 
free path. 

4. Summary and discussion 

We have calculated the time-dependent transport of UHECR in the Galaxy, assuming 
it can be described as isotropic diffusion. We have in particular investigated the question, 
whether or not Galactic long GRB can provide the observed particles up to the ankle, which 
is a natural transition point to extragalactic particles. In an earlier paper of the same series 
(Paper II) it was demonstrated that if extragalactic GRB are to produce UHECR above 
the ankle. Galactic GR should contribute a large flux of sub-ankle particles. Our present 
findings can be summarized as follows: 

• With the assumption of spherical symmetry, one overestimates the UHECR flux by a 
factor of a few, and likewise underestimates the anisotropy. 

• Intermittency becomes serious if the mean free path for scattering exceeds 100 pc, 
unless the source rate is much higher than 1 per Myr. On average. Galactic long GRB 
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Fig. 11. — Temporal variation of the anisotropy as in the bottom panel of Figure |9l but for 
a smaller mean free path a factor 5 smaller than assumed there, i.e. Amfp = 0.2 r^, at high 
energies. Again the average and the median at 2.4 EeV are indicated which are both close 
to 1%, approximately at the level of the Auger upper limits. 
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need to contribute only about 10^^ erg/s in accelerated particles to fully account for the 
observed particle flux at 10^^ eV, assuming a Bohmian mean free path at this energy. 

UHECR from Galactic long GRB can meet the observational limits on anisotropy only 
if the mean free path for scattering is sufficiently small. Contributing the observed sub- 
ankle particles (at 10^^ eV) requires Bohmian diffusion if the UHECR are as heavy as 
carbon. A light composition such as protons or helium requires sub-Bohmian diffusion, 
which is a highly unlikely situation for isotropic diffusion. We have not investigated 
the effects of a Galactic guiding field that may modify the probability of escape from 
the Galactic disk. 

Auger data suggest that at 10^* eV the composition is indeed light, thus posing a 
problem for the notion that Galactic GRB (or any other source class with similar pop- 
ulation statistic) produce the observed UHECR up to the ankle. (This measurement is 
not undisputed, though, for the KASCADE-Grande collaboration has just published 
their an alysis results which seem to fav or a relatively heavy composition up to nearly 
10^^ eV JArteaga- Velazquez et al.lboioh .) 



The UHECR composition is a very critical constraint, but its measurement is subject 
to considerable systematic uncertainties arising from its dependence on models for 
the development of air showers. It is imperative that measures be taken to better 
understand the air-shower physics near 10^^ eV. 

Much of the UHECR anisotropy arises from the expected location of long GRB in the 
inner Galaxy. Observations of GRB host galaxies suggest that region s of low metal 



licity and high star fo rmation may be the preferred sites of long GRB (ILevesque et al. 



2010l : iLevesqud l201ll ). which may skew the galactocentric distribution of long GRB 
toward the outer Galaxy. As there is no power problem with Galactic GRB, it may 
be worthwhile to also consider short GRB. They provide supposedly less power as 
a population, but they may have a very e xtended spatial distribution in the Galaxy, 



thus reducing the anisotropy (jBergerll201Cll ). The anisotropy arising from intermittency 



would remain in both cases and is the subject of a forthcoming publication. 

Our study can be applied with little change to the case of an UHECR origin in su- 
pernova remnants, assuming very efficient magnetic-field amplification can increase 
their ability to acc elerate particles to energies significantly higher than 1 PeV (e.g. 



Bell Sz Lucekl 120011 ). The spatial distribution in the Galaxy of long GRB and SNR 
can be expected to be similar, and therefore the average anisotropy is the same for 
both long GRB and SNR. We have verified that intermittency would be negligible in 
case of a SNR origin (unless of course the UHECR come only from a rare subclass of 
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supernovae) and therefore problems in simultaneously matching the composition and 
the anisotropy limits of UHECR are inescapable for SNR. 
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